Transport properties of 2D graphene containing structural defects 
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We propose an extensive report on the simulation of electronic transport in 2D graphene in pres- 
ence of structural defects. Amongst the large variety of such defects in sp^ carbon-based materials, 
we focus on the Stone- Wales defect and on two divacancy-type reconstructed defects. First, based 
on ab initio calculations, a tight-binding model is derived to describe the electronic structure of 
these defects. Then, semiclassical transport properties including the elastic mean free paths, mo- 
bilities and conductivities are computed using an order-N real-space Kubo- Greenwood method. A 
plateau of minimum conductivity {cr^c^ — ^e"^ /Tvh) is progressively observed as the density of defects 
increases. This saturation of the decay of conductivity to cr^*^ is associated with defect-dependent 
resonant energies. Finally, localization phenomena are captured beyond the semiclassical regime. 
An Anderson transition is predicted with localization lengths of the order of tens of nanometers for 
defect densities around 1%. 

PACS numbers: 73.23.-b, 72.15.Rn, 72.80.Vp, 71. 23. An 



INTRODUCTION 



Most exceptional properties of graphene are closely re- 
lated to its crystalline structure and its reduced dimen- 
sionality. The combination of a honeycomb lattice and a 
real two dimensional material is at the origin of the lin- 
ear dispersion and pseudospin symmetry of graphene low 
energy carriers. This gives rise to a wealth of interesting 
properties such massless Dirac fermion^, reduced back 
scattering, Berry phasei^ and weak antilocalization^. 
Besides, some transport characteristics of graphene, such 
as high carrier mobilities^ and long electronic coherence 
lengths^, are highly interesting from a technology point 
of view. Graphene could also match many request of the 
process and design of nanoelectronic devices. Indeed, be- 
cause of its 2D geometry, strategies of top-down fabrica- 
tion approach such as photolithography, intensively used 
in the mass production industry, could be more easily en- 
gineered and applied to design graphene-based nanode- 
vices architectures. Mainly for these reasons, graphene 
is expected to become a material of choice for future na- 
noelectronic. 

Yet, the lack of electronic (or transport) gap in pristine 
graphene is an issue that has to be overcome for achieving 
standard electronic devices such as field effect transistors 
(FET). As it has been the case for silicon, the tunning 
of graphene electronics is an essential step towards appli- 
cations. For this reason, controlled defect engineering in 
sp'^ carbon-based materials has become a topic of great 
excitement^. Indeed, the electronic (and transport) prop- 
erties of carbon nanotube^i and graphene-based materi- 
als^Tii^ can be considerably enriched by chemical modifi- 



cations, including substitution and molecular dopingiiii^ 
as well as functionalization. Another approach to tune 
graphene's properties is provided by ion or electron beam 
irradiation which introduces structural defects {e.g. va- 
cancies) in sp'^ carbon-based nanostructures. As an ex- 
ample, convincing room-temperature signatures of an 
Anderson regime in Ar^ irradiated carbon nanotubes 
have been reportedi^J^. In contrast, the conductivity 
of irradiated two-dimensional graphene saturates at the 
Dirac point above e'^/h even down to cryogenic tem- 
peratures^, suggesting a stronger robustness of defec- 
tive graphene. Nonetheless, many symmetries are broken 
when lattice disorder, such as structural defects, is intro- 
duced into a perfect honeycomb latticei^. Accordingly, 
structural defects in such materials should strongly affect 
the transport properties. We will show that a strong (An- 
derson) localization regime could be observed but, only 
in regions of energy corresponding to specific resonances 
directly associated to a certain type of defect and which 
not necessarily coincide with the Dirac point. Another 
reason why irradiated two-dimensional graphene seems to 
exhibit strong robustness compared to irradiated carbon 
nanotubes for an equivalent amount of structural defects, 
is the dimensionality. Indeed, the scaling theory of local- 
ization predicts a different behavior of the localization 
phenomena in one and two dimensional materials^^. 
The paper is elaborated as follow: In section [III we de- 
scribe the geometry and the electronic properties of struc- 
tural defects with ab initio and tight-binding techniques. 
Then, in section IIIIl we present the real-space Kubo- 
Greenwood methodology and discuss the transport prop- 
erties of large graphene planes with a realistic structural 
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FIG. 1. (color online). Schematic of the three structural defects: (a) Stone- Wales, (b) 585 and (c) 555-777 divacancies. 
Symmetry axis are drawn in red dashed lines. 



disorder. We also perform an in-depth analysis of the 
localization phenomena obtained within the present ap- 
proach. Finally, in section HVl we conclude on our results 
and discuss about the possibility to validate our theoret- 
ical predictions with experiments. 



II. ELECTRONIC PROPERTIES OF 
STRUCTURAL DEFECTS 

A. Geometry of structural defects 

Structural defects exist in various forms in 
graphenei^— and more generally in sp^ carbon- 
based materials. In this paper, three types of structural 
defects are considered, the Stone-Wales (SW) defect 
and two possible reconstructions of the divacancy. 
These three defects are non-magnetic, in contrast with 
monovacancies for instance which contain an unbounded 
carbon atom. Consequently, the issue of induced 
magnetization and spin interaction between defects does 
not need to be taken into account. The SW defect is 
a well known and common defect in sp^ carbon-based 
materials^^ which consists in a 90 degree rotation of a 
carbon-carbon bond. This topological transformation 
yields to the formation of two heptagons connected 
with two pentagons (FiglUa). Ma et al. have reported 
in a theoretical study^^ that in graphene a SW defect 
could produce a slight out-of-plane deformation. In this 
situation, the carbon bonds in the close neighboring 
of the SW defect are no more fully sp^ hybridized but 
become partially sp^-Yike hybridized. Considering this 
theoretical prediction would imply that the new orbitals 
hybridization could induce spin-orbit coupling as it has 
been predicted for the case of ripples^^. However, such 
out-of-plane deformations have never been obtained in 
our calculations. To ensure that the in-plane geometry 
was not a metastable state related to a local minimum 
of total energy, additional calculations have been per- 
formed with a starting geometry containing out-of-plane 
atoms. After structural relaxation, the out-of-plane 
atoms always go back to the plane. Moreover, to our 
knowledge no experimental STM measurements have 
confirmed this behavior predicted in ReS^. Therefore, 



only the in-plane geometry has been considered in this 
article. 

Vacancies are missing carbon atoms in the honeycomb 
lattice. They can not be considered as a reversible 
geometrical modification of the ideal graphene plane 
and thus strictly speaking, they do not belong to the 
class of topological defects. Vacancies can be created by 
irradiating the graphene plane with ions such as Ar^ for 
instance. Single vacancies (or monovacancies) migrate 
easily in the graphene plane and are stabilized when 
they recombine with another one to give a divacancy 
defectf^ i^^i^^ . Two kinds of divacancy defects are ex- 
plored here. In the first case, the reconstruction yields 
to the formation of 2 pentagons and 1 octagon (so-called 
585, Figdlb), while in the second case it yields to the 
formation of 3 pentagons and 3 heptagons (so-called 
555-777, Fig (He). According to our ab initio calculations 
the formation energy of the 555-777 divacancy is smaller 
than the one of the 585 divacancy by about 0.9eV. This 
stabilization of the 555-777 divacancy with regards to 
the 585 divacancy in graphene, which contrasts to the 
case of carbon nanotubes, is also reported by G.-D. Lee 
et al?^ . A third kind of divacancy, not studied here, has 
also been reported^^^^''. This third type of divacancy 
has a larger extension and is thus more complicated to 
parametrize within a tight-binding model. Actually, it 
involves 9 carbon rings including 4 pentagons, 1 hexagon 
and 4 heptagons (so-called 5555-6-7777). To conclude 
on the geometry analysis of these defects, one notes 
that the 585 divacancy as well as the SW defect possess 
a D2h symmetry since two orthogonal symmetry axis 
can be defined, whereas the 555-777 divacancy possess 
a D^h symmetry (see Fig{T]). The presence of these 
three structural defects has already been experimentally 
reported in graphene by means of STM experiment^i^S 
or TEM images2i, and their influence on the transport 
properties deserves in-depth inspection. 



B. Ab initio simulations 

In order to study these three defects, ab initio simu- 
lations are first employed. The calculations have been 
performed using the SIESTA package^^ in the general- 
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FIG. 2. (color online). Electronic band structures computed using the SIESTA package along high- symmetry lines in the 
Brillouin zone (left panel) for 7 x 7 supercell containing, (a) one SW defect, (b) one 585 divacancy, (c) one 555-777 divacancy. 
The Fermi energy is set to zero and the Dirac energy is indicated with a horizontal blue dashed line. 



ized gradient approximation (GGA) for the exchange- 
correlation functional in the Ceperley- Alder form^^, as 
parametrized by Perdew and Zunger^^. Troullier-Martins 
pseudopotentials are used to account for the core elec- 
trons^. The valence electron wave functions are ex- 
panded in a double- C polarized basis set of finite-range 
numerical pseudoatomic orbitals^^. 

The ab initio calculations have two objectives. On one 
hand, it allows to obtain the detailed geometry of the 
defect, and on the other hand, it allows to elaborate a 
tight-binding model for an isolated defect. In both cases, 
a convergence study with regards to the size of the cell 
used in calculations needs to be performed. The super- 
cell technique, which consists in placing the defect in a 
cell which is a multiple of the standard unitcell of the 
pristine system, is therefore employed. 
First, a structural relaxation of the supercell containing 
one defect is performed for different supercell sizes, i.e. 
4x4, 5x5,6x6 and 7x7(1x1 being the conventional 
unitcell). To ensure that the defect does not interact 
with its repeated images, the potential associated to the 
defect must tend to zero at the edges of the supercell. 
The 7x7 supercells are found to be large enough for 
each defect to obtain such a convergence. Final positions 
of the atoms in the vicinity of the defects are then ex- 
tracted and will be used to construct the tight-binding 
models. The 7x7 supercells after structural relaxation 
of each defect are displayed in Figdl The SW defect does 
not induce strong distortions of the neighboring hexagon 
rings, whereas the 585 divacancy exhibits a more pro- 
nounced impact. Indeed in the latter case, a distortion 
of the hexagons rings along the vertical axis (aligned with 
octagon) is observed. A slight distortion is also observed 
for the 555-777 divacancy but smaller than for the 585 
one. 

In FiglSl the electronic properties of the defects are ex- 
amined through the study of their respective band struc- 
tures. The electronic bands are plotted along an ex- 
tended path connecting high symmetry points in the Bril- 
louin zone. The path used here (see left panel of Figl2j), 
which is longer than the usual K-F-M-K path used for 
graphene, has been carefully chosen to evidence all use- 
ful details. Indeed, as the structural defects break the 



symmetry of the pristine graphene lattice, the hexagonal 
symmetry of the Brillouin zone is also lifted. As a conse- 
quence, paths in the Brillouin zone that were equivalent 
in the pristine case can become inequivalent when such 
a type of defect is introduced. By definition, if there are 
no other symmetries apart the time-reversal symmetry, 
half of the Brillouin zone needs to be examined which im- 
plies at least 10 branches (full and dotted red lines in left 
panel of Figl2j). Here, a minimal continuous path com- 
posed of 6 branches have been chosen in order to depict 
all important details for each defect. This extended path 
is labeled a-b-c-d-e-f in the left panel of Fig 121 
Then, for sake of comparison between band structures 
of the three defects, the Fermi energy has been set to 
zero and the energy associated with the Dirac point has 
been indicated by a horizontal blue dashed line. The 
Dirac cone band crossing occurs at an energy E = 0.7eV 
above the Fermi energy in case of the 555-777 divacancy 
(FigllJc), thus suggesting that this defect has a p-type 
doping character since the Fermi energy has been shifted 
below the Dirac point energy, whereas both energies were 
aligned in the pristine case. Obviously, the magnitude of 
the shift of the Fermi energy is related to the concentra- 
tion of defects (ud). One defect in a 7 x 7 supercell cor- 
responds to Ud ^ 3.9 lO^^cm"^ or equivalently rid ^ 1% 
which is the maximum concentration that will be consid- 
ered in section Hn] for transport properties. For the 585 
divacancy (Figl2lb), the Dirac point energy is located at 
E = 0.6eV above the Fermi energy, also implying a p- 
type doping character although slightly weaker than for 
the 555-777 divacancy. Finally the Dirac point energy 
and the Fermi energy are both aligned in case of a SW 
defect (Fig 121 a) which indicates a no doping effect. 
As mentioned in the previous paragraph, the 555-777 di- 
vacancy possesses a Dsh symmetry which is the same 
symmetry group as a single Dirac cone (at K or at K ). 
Indeed, even if very close to the Dirac point energy the 
Dirac cone is isotropic, at higher energies it has been 
demonstrated theoreticall y^^i^^ and by ARPES measure- 
ments^^ that a trigonal warping emerges. This effect is 
a signature of the Dsh symmetry of a single Dirac cone 
related to the existence of two inequivalent triangular 
sublattices. Since the 555-777 divacancy has the same 
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symmetry, the band structure of the defected system pre- 
serves the initial symmetry present in pristine graphene. 
Therefore, the usual electronic path composed of only 
three branches should be sufficient. It is actually easy 
to see for instance that in FiglSlc, paths r-K2 and F-K^ 
are equivalent as well as paths K2-K2 and K^-Ki. On 
the contrary, the 585 divacancy and the SW defect pos- 
sess a D2h symmetry which breaks the symmetry of the 
Dirac cone and as a consequence the latter is shifted. 
For the 585 divacancy, the first Dirac cone band crossing 
corresponding to K valley is shifted along the direction 
K2 ^ F, whereas the second Dirac cone band crossing 
corresponding to K valley is shifted in the opposite di- 
rection (K3 F). The latter is not visible in Figl2jb 
since the path K3-F is not represented (but is the sym- 
metric of K2-F). For the SW defect the shift is reversed, 
i.e. the first Dirac cone band crossing corresponding to 
K valley is shifted along the direction F ^ K2 (or equiv- 
alent ly along Ki K^), whereas the second Dirac cone 
band crossing corresponding to K valley is shifted along 
the opposite direction (K^ Ki). For both defects, 
paths F-K2 and F-K^ are now inequivalent as well as 
paths K^-K2 and K^-Ki. This shift of the Dirac cone in- 
duced by symmetry breaking has already been reported 
in case of uniaxial strain^, but to our knowledge not in 
case of point defects. 

Some features of the transport properties could already 
be anticipated from these electronic band structures anal- 
ysis. For instance, rather flat bands close to the Dirac 
point energy are observed. These flat bands have to be 
related to the presence of the defect, and can be seen 
as resonance energies associated with electrons localized 
around the defect. Such localized states are in general 
responsible of reduced transport properties. For the SW 
defect, a flat band is observed around E = 0.5eV above 
the Dirac point energy. A second flat band is located 
at ^ = — 1.75eV but is embedded in many dispersive 
bands. For the 585 divacancy, two regions of flat bands 
are observed around OeV and 0.25eV corresponding re- 
spectively to energies E ^ — 0.6eV and E ^ — 0.35eV 
below the Dirac point energy. Another flat band is also 
observed around — leV but again many other dispersive 
bands stands at this energy. Finally, for the 555-777 di- 
vacancy, much more flat band regions are observed but 
far away from the Dirac point energy. Indeed, flat bands 
are observed around 1.35eV, OeV and — leV which cor- 
responds respectively to energy E = 0.45eV above the 
Dirac point energy, E = — 0.7eV and E = — 1.7eV below 
the Dirac point energy. 

Finally, one concludes the analysis of band structures 
by investigating the case of charged systems. This is 
reported here as a validation of the rigid band model 
that will be used in the next section in order to discuss 
transport properties of defective graphene at different en- 
ergies. In contrast with the experimental configuration 
where the applied gate voltage effectively fill or deplete 
the active part of the device with electrons, our trans- 
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FIG. 3. (color online). Electronic band structures computed 
using the SIESTA package along high-symmetry lines F-M-K- 
F for 7x7 supercell containing, (a) one SW defect, (b) one 585 
divacancy, (c) one 555-777 divacancy, and with different ex- 
cess charge ({—3, —2, —1, +1, +2, +3} |e|). The Fermi energy 
is set to zero. 



port model relies on the rigid shift of the Fermi energy 
in order to describe out-of-equilibrium situations. In- 
deed, the tight-binding and transport models presented 
in the next paragraphs are essentially a model of the 
equilibrium case and can not account for the detailed 
impact of the excess of charge. In FiglH the validity 
of this rigid band approximation, which neglects these 
charge effects, is verified by looking at the impact of net 
charges on the electronic band structures. Calculations 
have been performed for systems with excess charge of 
{— 3, — 2, — 1, +1, +2, +3} |e| where |e| is the elementary 
charge. Note that, computationally, the added charges 
are compensated by a background jellium. For a 7 x 7 su- 
percell, a charge of 3 |e| corresponds to an carrier density 
of n ~ 1.2 lO^^cm"^. Figl3] shows that the overall band 
structure is not strongly modified by the added charges 
up to ±3|e| especially in the region [— l,+l]eV around 
the Fermi energy. For higher energies, the approxima- 
tion of rigid defect model is more questionable. Also, it 
is interesting to note that, when the conventional path 
composed of 3 branches (F-K-M-F) in the Brillouin zone 
is used to plot the band structure (see Figl3j), a fictitious 
energy gap is observed. This pinpoints the importance 
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of the symmetry breaking of the Brihouin zone and the 
need of such an extended path as used in Figl2l 



C. Tight-binding models 

The tight-binding (TB) models for pristine graphene 
and for the three defects are derived from the previous 
ab initio calculations. A common way to obtain the TB 
parameters is to choose a set of points E{k) in the ab 
initio (or experimental) band structure and then to use 
them as constraints in a fit procedure (based on gene al- 
gorithm for instance). The TB parameters are adjusted 
by the fit to reproduce the full band structure as well as 
possible. Here, a different strategy is applied. Since the 
SIESTA Hamiltonian is expressed in a localized orbitals 
basis set, the TB parameters can be directly extracted by 
performing successive operations on this SIESTA Hamil- 
tonian. In particular, the basis set has to be reduced to 
a single orbitals orthogonal basis. The details of the 
technique will be presented in a separated paper^^. 
Most TB studies use a nearest neighbors tt-tt* TB 
model to describe graphene electronic properties. How- 
ever, in the present work, a two centers 3^^ nearest neigh- 
bors TT-TT* orthogonal TB model has been chosen. In- 
deed, already for pristine graphene, a 3^^ nearest neigh- 
bors 7T-7T* model gives much more reasonable results since 
it allows in particular to recover the existing asymme- 
try between valence (tt) and conduction (tt*) bands. In 
contrast, a nearest neighbors tt-tt* model produces 
a totally symmetric band structure as shown in FigHJ 
The comparison between the ab initio and the 3^^ near- 
est neighbors model band structures and total densities 
of states (DOS) is very satisfactory. However, the K-M 
branch in the conduction band side remains not exactly 
reproduced by this 3^^ nearest neighbors model. The TB 
parameters of this model are composed of a single on-site 
term Sp^ and three hopping terms 7q, 7q and 7q corre- 
sponding respectively to l^^t 3rd nearest neigh- 
bors. The pristine graphene Hamiltonian then reads as 

{{4>^\7o\4>J) + {4>^\lo\<Pk) + {<P^\l|\<P^)) (1) 



DOSfeV am ) 



E 

i,<j,k,l> 



with Sp^ = 0.59745eV and 7^ = -3.09330eV, 7^ = 
0.19915eV, 7^ = -0.16214eV. The sum on index i run 
over all carbon orbitals. The sums over j, /c, / in- 
dexes run over all p^ orbitals corresponding respectively 
to l^st 2nd 3rd Clearest neighbors of the i^^ p^ orbital. 

To obtain the local TB parameters corresponding to 
the defect potential, the same extraction technique as for 
the pristine graphene have been used. In the defect po- 
tential, only on-site modifications have been considered, 
but the new arrangement of neighbors for carbon atoms 
in the core of the defects have also been carefully taken 
into account. In case of SW defect for instance, the rota- 
tion of the carbon-carbon bond yields to a modification 




FIG. 4. (color online). Electronic band structures and corre- 
sponding density of states (DOS) computed using the SIESTA 
package with a double-^" polarized (DZP) basis set (red lines) 
along K-F-M-K path for 1 x 1 supercell (unitcell). The TB 
band structures for a 1^^* nearest neighbors model (Inn, blue 
lines with open square symbols) and for a 3^^ nearest neigh- 
bors model (3nn, black lines with filled circle symbols) are 
also plotted. The Fermi energy is set to zero. 



of first, second and third nearest neighbors for carbon 
atoms in the vicinity of the rotated bond. To check the 
validity of the TB parametrization of the defects, the ab 
initio and the TB band structures have been compared 
for a 7 X 7 supercell containing one defect. In Figl5l the 
TB band structure (black dotted hues) is superimposed 
with the ab initio band structure (red full lines) already 
reported in Figl2l A good agreement is obtained espe- 
cially for the valence bands. The conduction band side 
seems to be less accurate but this is uniquely due to the 
inability of pristine graphene TB model to reproduce con- 
duction band along K-M branch, as already illustrated in 
FiglH As a consequence, some of the conduction bands 
are shifted to higher energies. A comparison of ab initio 
and TB DOS is also provided in the right panels of FigEl 



III. TRANSPORT PROPERTIES 

In this section, the transport properties of large defec- 
tive graphene planes is investigated. Although the elec- 
tronic structure study performed in the prior section al- 
ready allows to draw some preliminary conclusions on the 
impact of structural defects on transport, the randomness 
character of a reahstic disorder at the mesoscopic scale 
has to be taken into account. Therefore, based on the 
developed TB models for isolated defects, large graphene 
planes containing randomly distributed defects of vari- 
ous nature and density are built (Figl6j). It is important 
to pay attention to the various possible orientations of 
the defects. As explained previously, the structural de- 
fects possess specific symmetry groups smaller than the 
one of the honeycomb lattice. Thereby, for the SW and 
the 585 divacancy there exist three possible orientations, 
whereas there are only two for the 555-777 divacancy. If 
a given orientation is overrepresented, anisotropic trans- 



6 




K', r K', K, M,, r 4 8 12 16 



FIG. 5. (color online). Ab initio electronic band structure 
(red lines) computed using SIESTA with a double-C polarized 
(DZP) basis set along high- symmetry lines in the Brillouin 
zone (see left panel of Fig|2]) for 7 x 7 supercell containing, 
(a) one SW defect, (b) one 585 divacancy, (c) one 555-777 
divacancy. The TB band structures for a 3^^ nearest neigh- 
bors model (3nn, black lines with filled circle symbols) is also 
plotted. The Fermi energy is set to zero and the energy of the 
Dirac point is indicated with a horizontal blue dashed line. In 
right panels, the corresponding DOS are plotted. The pristine 
graphene DOS for a 3^^^ nearest neighbors model is superim- 
posed for sake of comparison (green thin lines). 



port properties would be observed. In our disorder model 
the defects are randomly oriented to avoid this effect. 
The case of graphene planes containing a mixture of de- 
fects is also studied in order to be as close as possible 
to the experimental situation. In a previous paper, the 
case of graphene planes containing only one type of de- 
fect but with different densities has been investigated^. 
For a mixture of several types of defects, the results can 
be almost derived by linear superposition of transport 
properties obtained for one specific type of defects as it 
will be demonstrated in the next paragraph. 




FIG. 6. (color online). Example of a small piece of defected 
graphene plane generated using the TB models developed for 
the three defects. Both position and orientation of the defects 
are chosen randomly. 



A. Density of states 

Before achieving an in-depth analysis of the transport 
properties, it is instructive to examine the density of 
states of randomly disordered graphene planes. Indeed, 
a first guidance to understand the transport features can 
be obtained by investigating the band structures and the 
DOS of the system as it has been performed in para- 
graph [nB] for defects in small supercells. Here, the DOS 
of these much larger graphene planes with randomly dis- 
tributed defects reveal the salient features that persist 
after taking into account the randomness character of 
the disorder. In Figd the total DOS of large graphene 
planes containing n^^ = 1% of SW, 585 and 555-777 di- 
vacancies, computed using the recursion methoc^, is 
compared with the total DOS previously obtained for a 
7x7 supercell containing one defect {ud ^ 1%). A first 
observation is that the DOS of random disordered sys- 
tems is much smoother than the one corresponding to a 
supercell^^. In the random case, most of the peaks have 
disappeared except the ones close to the Dirac point (In 
FiglTl the Dirac point has been set to zero.). The broad- 
ening due to the distribution disorder is more efficient 
in energy regions containing a lot of bands. Close to the 
Dirac point, the quantity of bands is less dense preserving 
the defect-induced resonances. Secondly, the position of 
resonance energy peaks corroborate the previous analysis 
carried out on supercell band structures. The presence 
of such resonant states induced by local defects is a well 
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FIG. 7. (color online). TB DOS for a 7 x 7 supercell (dashed lines) containing, (a) one SW defect, (b) one 585 divacancy, (c) 
one 555-777 divacancy, i.e. Ud ^ 1%, and for a large plane (thick lines) with Ud = 1% of, (a) SW defects, (b) 585 divacancies, 
(c) 555-777 divacancies, distributed and oriented randomly. 



known issueS^^^"—. The DOS of randomly disordered 
graphene confirm that the electron transport in an energy 
region around E = 0.35 eV is expected to be degraded 
for SW defects, whereas hole transport should be altered 
around E = —0.35 eV for 585 divacancies, and finally 
that 555-777 divacancies exhibit several resonance ener- 
gies around E = 0.6,-0.8,-2.1 eV which should also 
lead to reduced transport performances. 
The case of a mixture of defects is now explored by con- 
sidering graphene planes containing, half SW half 585 
divacancy (SW/585), then half SW half 555-777 diva- 
cancy (SW/555-777), and finally half 585 half 555-777 
divacancies (585/555-777). The corresponding DOS of 
these systems for n^^ = 1% of defects in total are plot- 
ted in FiglSl These DOS of graphene planes containing 
a mixture of two types of defects are compared also with 
the DOS of graphene planes containing a single type of 
defect separately (i.e. Ud = 0.5%). It is obvious that the 
features observed in the DOS of the mixed systems are 
roughly the sum of the individual features of each defect 
type. The particular case of SW/585 is interesting in 
the sense that the resonance peaks of these two defects 
are almost symmetric with respect to the Dirac point. 
Moreover, since the resonance peaks are relatively close 
to the initial Dirac point, they tend to overlap at this 
special point and yield to an increase of the DOS at the 
Dirac point. In this special situation, there is no more a 
clear minimum of DOS associated with the Dirac point. 
By adding other types of defects, a large increase of the 
DOS at the Dirac point can be foreseen. This is actu- 
ally what is observed for highly defective or amorphous 
graphene membranes^^^^. But this DOS increase comes 
from resonance states mainly localized around the defects 
which will therefore not participate to the transport of 
charge carriers but will rather degrade it. 



B. Kubo- Greenwood methodology 

The transport properties of large graphene planes con- 
taining structural defects are calculated using an efficient 
real-space order-N Kubo- Greenwood method^—. This 
very efficient technique gives a direct access to the main 
transport quantities in the semiclassical regime as well 
as in the quantum regime in which all multiple scat- 
tering events and interferences are retained. All quan- 
tities at energy E are extracted or derived from the 
wave packet dynamics. The latter is characterized by 
the time-dependent diffusivity D{E,t) = AR'^{E,t)/t 
where AR^ = AX^ + AY'^ is the mean quadratic spread- 
ing of wave packets. The mean quadratic spreading 
along a given direction is evaluated from the Hamiltonian 
and the corresponding position operator as AX'^{E^t) = 
Tr[(5(E - H)|X(t) - X(0)p]/Tr[(5(E - H)]. Tr is the trace 
over orbitals and Tr[(5(E - H)]/S = p{E) is the total 
DOS (per unit of surface). The two position operators 
X{t) and Y(t) are expressed in the Heisenberg represen- 
tation {X{t) = U^t)X{0)U{t)) and the time evolution 
operator U{t) = U^~q ex.p{iHAt/h), with At the cho- 
sen time step, is computed with a Chebyshev polyno- 
mial expansion method^ ~. Calculations are performed 
for several initial random phase wave packets, and for 
a total elapsed time t ~ 3.6ps split in three parts with 
different time steps, Ati ^ 1.32fs, At2 ^ 26.33fs and 
Ats ^ 52.66fs. The size of simulated graphene planes 
is X Ly = 300 X 250nm2 - 0.074/im2 (i.e. 2.8 x 10^ 
atoms), large enough to avoid finite size effects while peri- 
odic boundary conditions are applied to ensure continuity 
of the propagation of the wave packets. In the Kubo- 
Greenwood formalism, the different transport regimes 
can be inferred from the time dependence of the diffu- 
sivity coefficient D{E,t). The wave packet velocity v{E) 
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FIG. 8. (color online). TB DOS for a large graphene plane with Ud = 1% defects (thick solid lines) of, (a) SW/585, (b) 
SW/555-777, and (c) 585/555-777. These DOS are compared with DOS obtained with rid = 0.5% of each defects separately 
(dashed lines). The position of the Fermi energy is indicated by a vertical arrow. 



can be extracted for instance from the short time behav- 
ior of the diffusivity {ballistic regime)^ D{E^t) ~ v'^{E)t, 
while the elastic mean free path ie{E) is estimated 
from the maximum of the diffusivity {diffusive regime)^ 
Dmax{E) = 2v{E)£e{E). FoT two-dimensional systems, 
the Kubo- Greenwood conductivity equation is given by 

aiEF,TF) = - dE'^I^^^a{E\OK) (2) 

where is the Fermi energy and f{E^ Tp) is the Fermi- 
Dirac distribution function with Tp the associated Fermi- 
Dirac temperature. For Tp = OK, the Kubo-Greenwood 
conductivity is defined by 

a{EF,TF = OK) = ]e^p{EF) lim ^AR^{EF,t) (3) 

4 t^oo at 



classical conductivity 

(7sc{E,0K) = -e^p{E) lim D{E,t) (4) 

4 t^oo 

Gsc{E, OK) = \e^p{E)D^,^{E) (5) 

Finally, following the scaling theory*^— , the localiza- 
tion length i{E) can be determined from the semiclassi- 
cal transport length scales as, 

e(^) = V24(^)exp( -^--g'»^) ) (6) 

In the quantum regime, interferences in scattering paths 
survive and a decrease of D{t) is therefore expected {lo- 
calization regime). A sketch in Figl9] illustrates the typi- 
cal behavior of the diffusivity coefficient as a function of 
time, outlining the various possible transport regimes. 



In these two last equations, the Fermi level {Ep) has to 
be thought as an energy that can be tuned with a virtual 
electrostatic gate, so that the transport is accessible in a 
large energy window around the equilibrium position of 
the Fermi level. To avoid any confusion however, we set 
Ep = E in the following and reserve the notation Ep to 
the equilibrium Fermi energy at zero gate voltage (when 
the system is charge neutral). Note that we ignore the 
effects of screening by charge accumulation, which have 
been shown to be negligible in section HlBI 
In the semiclassical transport picture, the asymptotic be- 
havior of AR'^{E^t) at long time is linearly proportional 
to t {diffusive regime). Therefore, the derivative with re- 
spect to the time in Eql3] can be replaced by a simple 
division by the time. The diffusivity coefficient is thus 
introduced as D{t) = AR'^{t)/t in the Kubo-Greenwood 
conductivity formula to obtain the expression of the semi- 




max 



FIG. 9. (color online). Sketch of the typical behavior of D{t) 
showing the three possible transport regime including the lo- 
calization regime where the diffusivity decreases due to quan- 
tum interferences. In the semiclassical picture, D{t) stays in 
the diffusive regime [i.e. limt^oo D{t) = Dmax]. 
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FIG. 10. (color online). Mean free paths in graphene planes with, (a) SW/585 defects, (b) SW/555-777 defects, and (c) 
585/555-777 defects, for different defect concentrations ranging from Ud — 0.1% to 1.0%. 



C. Semiclassical transport - Diffusive regime 

The semiclassical transport regime excludes all quan- 
tum interferences effects and corresponding localization 
phenomena (weak and strong localization regimes) and 
simply captures scattering events within the diffusive 
regime, as for example derived within the self consistent 
Born approximation^!^. Therefore, to extract the semi- 
classical transport length scales, the Kubo-Greenwood 
formalism is restricted to the diffusive regime limit, i.e. 
the maximum of the diffusivity coefficient i^max as shown 
in Fig|9] (red line). In the following, the semiclassical 
transport properties of graphene planes containing vari- 
ous types and densities of structural defects will be de- 
scribed. First, the elastic mean free paths will be ex- 
amined, then the mobilities, and finally the semiclassical 
conductivities will be depicted. 

Mean Free Path 

In FigHOl the calculated mean free paths is presented 
for graphene planes in which SW/585 defects, SW/555- 
777 defects, or 585/555-777 defects have been incorpo- 
rated, with different defect concentrations ranging from 
rid = 0.1% to 1.0%. The mean free paths exhibit strong 
variations in energy, with dips associated to the defect 
resonance energies (bumps in the DOS). Note that from 
now, Fermi energy is always set to zero. The largest 
and fastest variation is obtained in an energy window 
corresponding to the 585 defect resonance energy, which 
appears just at the right side of the Fermi energy in case 
of SW/585 and 585/555-777 systems. The maximum dif- 
ference in mean free path is about one order of magni- 
tude. A rapid variation of the transport length scales 
with energy (or gate voltage) offer interesting perspec- 
tives for engineering mobility gaps^^. Actually, for the 
specific needs of logic devices which require an on and 
an off state, the case of 585/555-777 systems might be 
the most efficient one because the shape of the corre- 
sponding curve £e{E) is the closest to a stepwise func- 
tion. For Ud = 0.1%, the maximum calculated mean 
free paths for all different defective systems are in the 



range of [60, 200] nm whereas the minimal ones lie within 
[2, 10] nm. For Ud > 1%, the mean free paths are found to 
be smaller than lOnm for any energy in the window con- 
sidered here {i.e. [—3, +3]eV). Finally, regarding the vari- 
ation of the mean free path with defect densities, better 
ratios between minimum and maximum mean free paths 
are expected for defect concentrations Ud < 0.1%. Unfor- 
tunately, smaller defect concentrations are not accessible 
to our simulations due to computational limitations. 

Mobility 

In FigdU the mobilities as a function of carrier en- 
ergy (upper panels) and of carrier density (lower panels) 
computed with a Fermi-Dirac temperature of 300K are 
illustrated. The mobility (/i) is calculated using the con- 
ventional definition ii{E^Tf) = cFsc{E.,Tf) / en^E^Tp) 
where e is the elementary charge and n{E^TF) the charge 
carrier density (electrons or holes) computed by integrat- 
ing the DOS from the Dirac point to the energy E. Since 
/i is inversely proportional to the charge carrier density 
(n), there is an obvious mathematical divergence of /i 
when n tends to zero (close to the Dirac point) as ob- 
served in FigfTTl In each of the three cases (SW/585, 
SW/555-777, and 585/555-777), the p-doping character 
induced by divacancies is clearly illustrated by the con- 
tinuous shift of the mobility peak with respect to the 
Fermi level as Ud increases (upper panels). As for the 
mean free paths, some dips are also visible even if they 
seem not be as deeper because the logarithm scale used 
in FigiUlrun over five orders of magnitudes. The fastest 
variation is again obtained around the resonance energies 
associated to 585 defects and yields to a change of mobil- 
ity values between two and three orders of magnitudes. 
For the largest defect concentration {ud = 1%), the mo- 
bility can drop below 10cm^V~^s~^ but keeping a value 
of a few ten thousands of cm^V~^s~^ at the Dirac point. 
In the lower panels of FigHH the mobilities are plotted 
versus the electron and hole carrier densities measured 
with respect to the Dirac point. Therefore, the mobility 
peaks are all aligned. Variations with the carrier density 
are found to be smoother than with the carrier energy. 
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FIG. 11. (color online). Mobilities as a function of charge carrier energy (upper panels) and as a function of charge carrier 
density (lower panels) of graphene planes with, (a)(d) SW/585 defects, (b)(e) SW/555-777 defects, and (c)(f) 585/555-777 
defects, for different defect concentrations ranging from Ud — 0.1% to 1.0%. 




FIG. 12. (color online). Conductivities as a function of charge carrier energy (upper panels) and as a function of charge carrier 
density (lower panels) of graphene planes with, (a)(d) SW/585 defects, (b)(e) SW/555-777 defects, and (c)(f) 585/555-777 
defects, for different defect concentrations ranging from rid — 0.1% to 1.0%. 
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which is due to the non hnear dependence between n 
and E {n{E) ex E'^). Nevertheless, the variation of hole 
mobility close to the Dirac point is reasonably sharp in 
case of 585/555-777 defective systems. Far away enough 
from the Dirac point, the mobility generally tends to a 
constant value meaning that the conductivity is linearly 
dependent to the charge carrier in these energy regions. 

Semidassical Conductivity 

In Fig (121 the semidassical conductivities (age) are 
plotted as a function of carrier energy and density and for 
Tp = 300K. As for mean free paths and mobilities, the 
conductivities exhibit dips at resonance energies. The 
value of conductivity is minimal at those energies, which 
indicates that a minimal conductivity can be obtained 
for an energy different to the Dirac point. There is also 
a minimum value of conductivity associated to the Dirac 
point which is visible for Fermi-Dirac temperature equal 
to OK (although not shown here). This minimum is very 
sharp and thus easily broadened by the Fermi-Dirac tem- 
perature. Nevertheless, for SW/555-777 systems, this 
minimum associated to the Dirac point is still observable 
in FigUlb for Ud = 0.1% {E - 0.2eV). 
When the defect concentration increases, the conductiv- 
ity at the defect resonances reaches the semidassical limit 
of minimum conductivity cr^*^ = 2Go/7r = 4e^/7r/i, 
leading to the progressive formation of a plateau of con- 
ductivity (see ref^^ for a discussion on the minimum 
conductivity). Consequently, the conductivity curves 
present a lot of different behavior versus the energy. 
These different behaviors (constant, linear, quadratic,...) 
are even more obvious in the lower panels of FiglT2l where 
conductivities are plotted versus the carrier density. In- 
deed, for Ud = 0.1%, depending on the defective system 
and on the type of carriers, either a linear behavior can 
be observed (for electrons of SW/585 and in a smaller 
extend for electrons of SW/555-777), or a quadratic be- 
havior (for electrons and holes close to the Dirac point 
and more precisely around the defect resonances of all 
three type of defective systems), or again a sub-linear 
behavior (for holes of SW/585 and electrons of 585/555- 
777). For Ud = 1.0% of defects, the conductivity even 
displays an almost constant behavior over a large range 
of carrier density since the conductivity has reached the 
semidassical conductivity limit cr^*^. The fact that all 
these different behaviors for asc{n) are present in these 
defective graphene systems is an interesting observation 
for the current puzzling interpretation of <Tgc(n) ^^i^^ . 
The semidassical electronic transport in graphene con- 
taining various mixtures of structural defects has been 
investigated in this paragraph. The structural defects, 
which break locally the symmetry of the honeycomb lat- 
tice, have been found to yield to resonance scattering 
at specific energies. These resonant defect levels induce 
peaks in the DOS which increase with the defect density. 
These resonant states conduct systematically to a degra- 
dation of all transport properties examined (mean free 
paths, mobilities, and conductivities). This degradation 



occurring at specific energies, constitute a fingerprint re- 
lated to each type of defect separately even when several 
types of defects coexist in the system (provided that res- 
onant levels are not to close in energy to overlap, which 
could be the case with other defects not studied here). 
The 585 divacancy produces a rather sharp resonance 
compared to the two other defects. In particular, the case 
of 585/555-777 defective system presents an interesting 
profile regarding variations of transport properties with 
respect to carrier energy and carrier density. Finally, 
large plateaus of minimum conductivity are predicted for 
sufficiently high concentration of defect {ud > 0.5%). 
In the next section, the quantum regime of transport in 
such defective graphene systems will be discussed. In this 
transport regime, the quantum interferences in electron 
trajectory will be taken into account which yields to lo- 
calization phenomena and thus an additional degradation 
of transport properties such as the conductivity. 

D. Quantum transport - Localization regime 

In the localization transport regime, quantum interfer- 
ences can yield to an enhanced probability to be back- 
scattered. In a common picture, this effect can be viewed 
as constructive interferences between forward and back- 
ward electron scattering trajectories that form loops. 
This is possible only if the electron keeps its phase and 
thus, any source of phase decoherence such as electron- 
phonon scattering will destroy these interferences. This 
is the reason why such localization effects are usually ob- 
servable at low temperature only. However, in graphene 
the electron-phonon scattering have been found to be ex- 
tremely low^i^, allowing the localization phenomena to 
occur even at room temperature. Concerning the trans- 
port quantities, the quantum interferences induce a cor- 
rection which can be described theoretically trough a 
Cooperon teru^. For a two-dimensional system, this 
Cooperon correction yields for instance to a decrease of 
conductivity which scale logarithmically with the sam- 
ple size (weak localization regime). For even stronger 
disordered systems, the scaling of the conductivity ex- 
hibits an exponential decay (strong localization ergime), 
and the electronic states become strongly localized such 
that the system behaves as an insulator from a trans- 
port properties point of view. Several theoretical and 
experimental studies on graphene have reported such a 
(semi) metal-insulator transition through a Anderson lo- 
calization of electronic states which was cause either by 
hydrogen atom^Si— , epoxy oxygen atoms^^Ti^, vacan- 
cies and other structural defects^Sj^i^ — , or even noble 
metal cluster^. 

Diffusivity 

In our real-space Kubo-Greenwood calculations, this 
degradation of transport quantities caused by the quan- 
tum interferences can be directly observed through the 
time dependence of the diffusivity coefficient. The 
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FIG. 13. (color online). Normalized diffusivity as a function of time for four selected energies (upper panels), and 3D plots as 
a function of time and energy (lower panels), of graphene planes with, (a)(d) SW/585 defects, (b)(e) SW/555-777 defects, and 
(c)(f) 585/555-777 defects, for a defect concentration Ud = 1.0%. 



Figs|T3la-c illustrate this time dependence of normal- 
ized diffusion coefficients {D{t)/ Dmax) for four selected 
energies (including the Fermi energy), while FigslT3ld-f 
present a 3D plot of the same quantity for a large energy 
window of the spectrum. The four energies selected in the 
Figs (131 a-c are pointed with vertical arrows in FigsIT3ld- 
f. The decreasing of D{t) clearly reveals the emergence of 
localization effects. This decreasing is not homogeneous 
in the whole spectrum but stronger at the defect reso- 
nance energies. Actually, far away from these resonance 
energies, there is no more decreasing and the behavior 
of D{t) at long time is a constant which is typical of a 
semiclassical diffusive regime. 

Localization Length 

Since localization phenomena at some specific energies 
have been observed, the corresponding localization 
length (^) can be estimated using Eql6] (see FiglT4|). For 
the highest defect concentration (1.0%), the smallest 
localization lengths, obtained exactly at the defect 
resonant energies, are ^ = 10-20nm about. At these 
energies, the conductivity reaches the semiclassical 
lower limit {i.e. cr^*^), and therefore, Eql6] becomes 
^{E) = V2ie{E) eyip{2) = 10.5ie{E) with the smallest 
value obtained for being around Inm. Away from 
the resonance energies, the localization length increases 
very rapidly mostly driven by the conductivity which 
appears in the exponential function. This exponential 
variation with the conductivity is a common feature of 
two-dimensional systems which make the experimental 
observation of localization generally difficult even if the 
mean free path is very small. 



Extension of the real-space Kubo- Greenwood 
methodology 

In the section llll C\ the Kubo-Greenwood methodology 
has been used to describe semiclassical transport prop- 
erties using quantities extracted only from the diffusive 
regime (maximum of D{t)). In the previous paragraph, 
even the localization lengths have been extracted from 
semiclassical quantities using the scaling theory of local- 
ization (Eql6]). The reason why the Kubo-Greenwood 
formalism is usually restricted to the diffusive regime, is 
because this is a hulk formalism in the sense that the 
framework is devoted to the description of the intrin- 
sic properties of a material. For instance, within the 
Kubo-Greenwood formalism, the calculated semiclassi- 
cal conductivity is an intensive quantity contrarily to 
the conductance which is extensive and depends on the 
system geometry. Hence, the Kubo-Greenwood conduc- 
tivity formula is usually employed in its thermodynamic 
limit t ^ oo (EqsISJH)). This infinite time limit imposes 
to restrict the analysis of D(t) to its asymptotic value 
which can be either a constant (purely diffusive regime, 
o' = Csc) or zero (completely localized regime, a = 0). 
On the contrary, the Landauer-Biittiker transport ap- 
proach, is better suited to describe extensive transport 
properties such as the conductance {G{L)). In this last 
formalism, the presence of leads introduces a finite length 
L of the system which allows scaling analysis describing 
therefore all intermediate regimes. Even if there is a pri- 
ori no clear definition of length in the Kubo-Greenwood 
approach since the system is periodic (bulk), the present 
real-space implementation of this transport formalism 
provides a complete description of the dynamics of prop- 
agation, going from the ballistic regime at very short 
time simulations to quasi- diffusive and diffusive regimes 




as well as the weak and strong localization regimes at 
longer time. In fact, in these real-space simulations, the 
computed spreading of wave packets can be seen as an 
average time-dependent distance probed by the propa- 
gating charge carriers. To allow a scaling analysis within 
the real-space Kubo-Greenwood approach, the introduc- 
tion of a characteristic length scale L(t) defined as the 
diameter of a fictitious circle formed by the spreading of 
wave packets is therefore proposed. 

L{E,t) = 2^AR^{E,t) (7) 

This new approach, which will give access to length 
dependent transport quantities and allow to describe 
more in details the localization phenomena, will be 
verified and justified in the following by comparing the 
results with the predictions of the scaling theory of 
localization. Before that, the definition of the diffusion 
coefficient D{t) has to be reexamined. 

As pointed out by Eql3l the general definition of the 
diffusion coefficient is given by D{t) = §^/^R^{t). In the 
diffusive regime, the quadratic spreading of wave pack- 
ets being linearly proportional to the time {AR'^{t) ex t), 
its first derivative reduces exactly to the constant value 
D{t) = AR'^(t)/t. For numerical convenience, it is more 
efficient to use the definition D{t) = AR^{t)/t. In- 
deed, the first derivative is more computationally de- 
manding and its numerical evaluation can produce large 
fiuctuations from which it is difficult to extract semi- 
classical transport properties. In the localization regime, 
since D{t) is no more a constant, the general definition 
D(t) = -^AR^it) should be used. However, although 
using D{t) = AR^(t)/t is in principle incorrect in the 
localization regime, it is nevertheless interesting to com- 
pare the two definitions. In Fig(T5l the diffusivity co- 
efficient is plotted using either D{t) = AR^{t)/t (solid 
lines) or D{t) = ^AR^{t) (dashed lines). For sake 
of comparison, both definitions are normalized with the 
same value of maximum of diffusivity (I^max) found for 
D(t) = AR^{t)/t. For an energy far from the resonance 
energies, e.g. E = —3eV for which no localization effects 



are observed, the two curves (solid and dashed lines) 
are almost superimposed. Both evaluations of D{t) at 
long time, deep in the diffusive regime, gives thus the 
same constant value. The only difference is a bump 
obtained for the diffusivity evaluated with the numeri- 
cal derivative. This bump appears at the transition be- 
tween the ballistic and the diffusive regime. The bump 
is followed by a small decrease of D{t) before the sat- 
uration is obtained. In principle, there is no reason to 
observe a decreasing of D{t) in a purely diffusive regime 
as shown in FiglQl This bump is thus just an artifact 
produced by small numerical errors in AR^{t) then en- 
hanced by its first derivative calculation^^. In the lo- 
calization regime, which is clearly predicted for energies 
E = -0.25, 0,0. 35eV, Fig US] shows that the diffusivity 
computed with D{t) = AR^(t)/t overestimates the one 
derived with D{t) = ^AR^(t). Nevertheless, the overes- 
timation in the localization regime being almost constant 
with increasing time, the general behavior is still correct 
when D{t) is evaluated with D{t) = AR^{t)/t. More- 
over, the curves obtained with that definition are much 
smoother. The definition D{t) = AR^{t)/t could there- 
fore be conserved even in the localization regime, keep- 
ing in mind that it slightly overestimates the true value. 
With such a clarification for the definition of diffusiv- 
ity in the localization regime, the scaling analysis of the 
conductivity in the quantum regime can be performed. 

By using Eq|71 transport quantities can be expressed as 
a function of length (L) instead of as a function of time 
{t). To check the validity of this presumed definition 
of L, the scaling of the conductivity can be compared 
to the predictions of the scaling theory. In particular, 
this theory predicts the length dependence of quantum 
corrections (Cooperon term) to the conductivity due to 
localization effects. For two-dimensional systems, these 
corrections are logarithmic^ and read as. 



(j{E,L)=(Jsc{E)- 

= (Jsc{E)- 
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FIG. 15. (color online). Normalized diffusivity D{t) in 
graphene planes with 1.0% of 585/555-777 defects, as a func- 
tion of time for five selected energies. Solid lines are ob- 
tained with Ai?^ (t) / 1 whereas dashed lines are obtained with 
■^AR^{t). Both expressions are divided by Dma^ found for 
AR'{t)/t. 
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localization effects induce a decreasing of quantum con- 
ductivity with increasing length. The semiclassical value 
asc (which is maximum of cr{L) = amax) is always found 
to be greater than cr^*^ = 2Go/7r, while quantum con- 
ductivity can decrease below such a value for sufficiently 
long lengths and for energies close enough to the defect 
resonance energies. Finally, Fig (16] exhibits a maximum 
computed length which is different for each energy. This 
is explained by the length scale (Eq|7j) which is energy 
dependent, thus for a fixed total elapsed time, the length 
probed by the spreading of wave packets can be different. 
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FIG. 17. (color online). Quantum conductivity in graphene 
planes with 0.25% of 585/555-777 defects, as a function of 
length and for three selected energies. Solid lines are obtained 
usms; D{t) = AR^(t)/t. The dotted lines are obtained usine: 
Eq|8]the and dashed-dotted lines using Eq|9]in which the new 
length scale /diff has been introduced. 
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FIG. 16. (color online). Quantum conductivity a{E,L) in 
graphene planes with 1.0% of 585/555-777 defects, as a func- 
tion of length for five selected energies. We use D{t) = 
AR^{t)/t (solid lines) and D{t) = ^AR^{t) (dashed lines) for 
computing cr(L). The horizontal dotted line gives the semi- 
classical limit a^J"^ = 2/71G0 = 4e^/7r/i. 



In Fig(T6l using Eqd the quantum conductivity as 
a function of length {cr{L)) is presented employing ei- 
ther D{t) = AR^{t)/t (solid lines) or D{t) = §iAR'^{t) 
(dashed lines). The horizontal dotted line denotes cr^^^. 
As for the diffusivity (FiglTSj) both definitions of D{t) 
lead to similar conductivity curves but the use of D{t) = 
AR^{t)/t produces smoother curves and a slight overes- 
timation in the localization regime. Localization effects 
are absent for E = — 3eV as evidenced by the saturation 
of (j{L) to its semiclassical asymptotic constant value cisc- 
For other energies (especially for E = — 0.25, 0, 0.35eV), 



In FiglTTl the validity of the definition of the length L 
(Eqd) is checked by scrutinizing if the quantum conduc- 
tivity decay (solid lines) follows a logarithmic behavior as 
predicted by the scaling theory (Eql8]). The dotted lines 
correspond to Eq|8l The dashed-dotted lines correspond 
to a modified version of Eql8]in which a new length scale 
(/diff) has been introduced replacing the mean free path, 

with /diff defined as the length corresponding to (Jmax- 
In the review of Lee and Ramakrishnan^^, the equation 
of weak localization correction in the scaling theory is 
derived using the mean free path as the characteristic 
length scale of the diffusive regime. According to this 
definition, when the length L is equal to the mean free 
path ^e, the conductivity is maximal and corresponds 
to the semiclassical conductivity {asc = cTmax), and for 
L > ie weak localization develops. In contrast, using 
Eq|71 the obtained semiclassical conductivity ( or amax) 
occurs at a length scale different from i^. For instance 
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in FigdTl for an energy E = OeV {E = Ef)^ the maxi- 
mum of conductivity occurs at L = /diff ^ 53nm, which 
is about ten times greater than the mean free path com- 
puted at this energy {i^ ~ 4nm). Therefore, the def- 
inition of length EqlTl might be questioned and argued 
to be misleading. But, considering the scaling in time 
for which there is no ambiguity, there is also a factor of 
ten between the calculated mean scattering time r and 
the time t for which the diffusivity coefficient reaches its 
maximum value. Actually, in a semiclassical picture, this 
is consistent with the fact that to reach completely the 
diffusive regime (saturation of D{t))^ the charge carriers 
needs to encounter several scattering events which corre- 
sponds thus to a total traveling distance of a few times 
the mean free path. The same idea also explains that 
the charge carriers need several scattering events before 
the weak localization corrections start to present a clear 
effect. Indeed, these corrections are due to constructive 
interferences between forward and backward scattering 
paths which form loop trajectories. The shortest loop 
trajectory is at least two times the mean free path in a 
semiclassical picture. The introduction of an additional 
length scale /diff is therefore meaningful for the interpre- 
tation of such quantum corrections. In Fig(T71 it is obvi- 
ous that the conductivity decay follows the logarithmic 
behavior predicted by Eq|8l However, using the mean 
free path as a characteristic length scale in the formula 
does not allow to fit cr(L). On the contrary, Eql9] allows 
to more correctly predict the scaling of the conductiv- 
ity and therefore to extrapolate to lengths L longer than 
the ones computed. As a consequence, the localization 
length ^ can be extracted with a better precision by re- 
placing the mean free path ^e by the diffusion length /diff- 
Indeed, ^ can be evaluated as the length L for which the 
conductivity is equal to zero or equivalently for which the 
weak localization corrections are equal to the semiclassi- 
cal conductivity. 

a{E, L = = = cT,e(i?) - ^ In f 

m = lME)e^v{"-^^) (10) 

Compared to EqlHl the localization lengths calculated 
using Eq(TO]are roughly ten times greater for weak local- 
ization regime and about 4-5 times greater for stronger 
localization regime. For instance, in Fig(T71 for an en- 
ergy E = 0.25eV which corresponds exactly to the res- 
onance of 585 divacancy and for which the localization 
effects are stronger, ^ ~ 25nm obtained with Eql6] and 
^ ~ lOOnm estimated with EqlTOl Note however that in 
case of strong localization effects such as ^ = 0.25eV in 
FigdTl and for higher defect densities such as 1% (not 
shown here), there is a deviance between the logarith- 
mic decay predicted by scaling theory for the weak lo- 
calization corrections and the actual computed quantum 
conductivities. As already reported recently ir^, in such 
strong localization case, the weak localization corrections 



are no more suitable and an exponential decay has to be 
used to describe the scaling of cr(L) which results in long 
tail corresponding to evanescent mode. Finally, in very 
strong localization regime, the localization length could 
alternatively be evaluated from the long time saturation 
of the quadratic spreading^ as 

£^{E) = lim 2y/AR'^{E,t) = lim L{E,t) (11) 

Unfortunately, in our simulations, a complete saturation 
of the quadratic spreading has never been observed for 
the total elapsed time computed. 



IV. CONCLUSION 

In conclusion, the electronic and transport properties 
of graphene planes containing structural defects such as, 
SW, 585 and 555-777 divacancies, has been investigated 
theoretically. First, ab initio techniques have been used 
to study the geometry and the electronic band struc- 
tures of single defects in a supercell. From these re- 
sults, relaxed atomic positions and electronic defect po- 
tentials have been employed to build accurate TB mod- 
els for each individual structural defects. Comparisons 
between ab initio and TB band structures and DOS 
have shown a good agreement. The elaborated TB mod- 
els have then been employed to compute DOS of large 
graphene planes containing randomly distributed struc- 
tural defects of various natures and concentrations. The 
comparison of these DOS with the DOS obtained for su- 
percell have put in evidence the smoothening effect of 
the randomness character in distribution and orientation 
of defect at a mesoscopic scale. Then, the computed 
DOS of graphene planes including a mixture of defects 
have illustrated that defect-induced resonant peaks can 
be assigned to individual defect separately demonstrat- 
ing that those resonances constitute a fingerprint of each 
structural defect. Afterwards, using a Kubo- Greenwood 
method, the transport properties have been examined. 
The computed semiclassical transport quantities, such 
as mean free path, mobility, and conductivity have ex- 
hibited systematic degradations around resonant defect 
energies. In particular, it has been predicted that a min- 
imum conductivity (cr^*^ = 4e^/7r/i) can be obtained at 
an energy different than Dirac point, and that plateau of 
minimum conductivity are formed when increasing the 
defect density. Finally, a detailed analysis of the diffu- 
sivity deep into the localization regime has allowed to ex- 
tend the use of the Kubo-Greenwood method to predict 
length dependence of the quantum conductivity (cr(L)). 
The latter exhibits a logarithmic decay with increasing 
lengths around defect resonance energies due to local- 
ization phenomena as expected by the scaling theory. 
An Anderson insulating behaviour is thus anticipated for 
graphene with 1% of structural defect and associated lo- 
calization lengths are expected to be in the order of a 
few tens nanometers. However, localization effects are 
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not observed in all the spectrum but only around defect 
resonant energies complicating therefore its experimental 
observation. Transport measurements would thus require 
an efficient electrostatic gate to probe larger part of the 
spectrum and to align Fermi energy with these defect res- 
onant energies. Additionally, transport characterization 
with different sizes of graphene samples exposed to the 
same amount of structural defects could be conducted 
in order to perform scaling analysis of the conductivity. 
Alternatively, magneto-transport experiments on a fixed 
size sample could also reveal localization phenomena in- 
duced by structural defects as well as low temperature 
transport measurements which should evidence a vari- 
able range hopping behaviour. 



J.-C.C. and A.L. acknowledge financial support from 
the F.R.S.-FNRS of Belgium. S.M.-M.D. acknowledges 
funding from EPSRC (Grant n.° EP/G055904/1), X.D. 
acknowledges financial support from the FRIA. Parts 
of this work are connected to the Belgian Program on 
Interuniversity Attraction Poles (PAI6), to the ARC 
Graphene Nano-electromechanics (n.° 11/16-037) spon- 
sored by the Communaute Francaise de Belgique, to 
the European Union through the ETSF e-I3 project 
(grant n.° 211956), and to the NANOSIM-GRAPHENE 
Project (n.° ANR-09-NANO-016-01). Computational re- 
sources have been provided by the CISM of the Universite 
catholique de Louvain (UCL). 



^ K.S Novoselov, Rev. Mod. Phys. 83, 837 (2011). 

^ Y. Zhang, Y.-W. Tan, H.L. Stormer, and P. Kim, Nature 

438, 201 (2005). 
^ F. Ortman, A. Crest i, G. Montambaux, and S. Roche, Eur. 

Phys. Lett. 94, 47006 (2011). 
^ K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. 

Hone, P. Kim, and H.L. Stormer, Solid State Comm. 146, 

351 (2008). 

^ F. Miao, S. Wijeratne, Y. Zhang, U.C. Coskun, W. Bao, 

and C.N. Lau, Sience 317, 1530 (2007). 
^ A.V. Krasheninnikov, and F. Banhart, Nature Materials 

6, 723 (2007). 

J.-C. Charlier, X. Blase, and S. Roche, Rev. Mod. Phys. 
79, 677 (2007). 

^ K. Suenaga, H. Wakabayashi, M. Koshino, Y. Sato, K. 

Urita, and S. lijima. Nature Nanotechnology 2, 358 (2007). 
^ M.T. Lusk, D.T. Wu, and L.D. Carr, Phys. Rev. B 81, 

155444 (2010). 

^° A. Cresti, N. Nemec, B. Biel, G. Niebler, F. Triozon, G. 
Cuniberti, and S. Roche, Nano Res. 1, 361 (2008). 

A. Lherbier, X. Blase, Y.M. Niquet, F. Triozon, and S. 
Roche, Phys. Rev. Lett. 101, 036808 (2008). 

^2 B. Wang, and S.T. Pantelides, Phys. Rev. B 83, 245403 
(2011). 

C. Gomez- Navarro, P.J. De Pablo, J. Gomez-Herrero, B. 
Biel, F.J. Garcia- Vidal, A. Rubio, and F. Flores, Nature 
Materials 4, 534 (2005). 

B. Biel, F.J. Garcia- Vidal, A. Rubio, and F. Flores, Phys. 
Rev. Lett. 95, 266801 (2005). 

J.-H. Chen, W.G. Cullen, C. Jang, M.S. Fuhrer, and E.D. 
Williams, Phys. Rev. Lett. 102, 236805 (2009). 
F. Libisch, S. Rotter, and J. Burgdorfer, 
arXivTl 102.3848 vl (2011). 
^'^ A. Lherbier, B. Biel, Y.M. Niquet, and S. Roche, Phys. 
Rev. Lett. 100, 036803 (2008). 

F. Banhart, J. Kotakoski, and A.V. Krasheninnikov, ACS 
NANO 5, 26 (2011). 

J. Kotakoski, A.V. Krasheninnikov, U. Kaiser, and J.C. 

Meyer, Phys. Rev. Lett. 106, 105505 (2011). 
^° E. Cockayne, G.M. Rutter, N.P. Guisinger, J.N. Grain, 

P.N. First, an d J. Stroscio, Phy s. Rev. B 83, 195425 (2011). 
2^ E. Cockayne. [arXiv:1106.6273V 2 (2011). 

A.J. Stone, and D.J. Wales, Chem. Phys. Lett. 128, 501 

(1986). 



2^ J. Ma, D. Alfe, A. Michaelides, and E. Wang, Phys. Rev. 

B 80, 033407 (2009). 
'^^ D. Huertas-Hernando, F. Guinea, and A. Brataas, Phys. 

Rev. B 74, 155426 (2006). 
2^ G.-D. Lee, C.Z. Wang, E. Yoon, N.-M. Hwang, D.-Y. Kim, 

and K.M. Ho, Phys. Rev. Lett. 95, 205501 (2005). 
2^ Y. Kim, J. Ihm, E. Yoon, and G.-D. Lee, Phys. Rev. B 84, 

075445 (2011). 

J.C. Meyer, C. Kisielowski, R. Erni, M.D. Rossell, M.F. 
Crommie, and A. Zettl, Nano Lett. 8, 3582 (2008). 
2^ M.M. Ugeda, L Brihuega, F. Hiebel, P. Mallet, J.-Y. 
Veuillen, J.M. Gomez- Rodriguez, and F. Yndurain, Phys. 
Rev. B 85, 121402(R) (2012). 

J.M. Soler, E. Artacho, J.D. Gale, A. Garcia, J. Junquera, 
P. Ordejon, D. Sanchez- Portal, Journal of Physics: Con- 
densed Matter 14, 2745 (2002). 

D. M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 
(1980). 

J.P Perdew, and A. Zunger, Phys. Rev. B 23, 5048 (1981). 
N. Troullier, and J.-L. Martins, Phys. Rev. B 43, 1993 
(1991). 

E. Artacho, D. Sanchez- Portal, P. Ordejon, A. Garcia, and 
J.M. Soler, Phys. Status Solidi. B 215, 809 (1999). 

2^ E. McCann, K. Kechedzhi, V.L Fal'ko, H. Suzuura, T. 
Ando, and B.L. Altshuler, Phys. Rev. Lett. 97, 146805 
(2006). 

M. Mucha-Kruczynski, O. Tsyplyatyev, A. Grishin, E. Mc- 
Cann, V.L Fal'ko, A. Bostwick, and E. Rotenberg, Phys. 
Rev. B 77, 195403 (2008). 

S.Y. Zhou, D.A. Siegel, A.V. Fedorov, and A. Lanzara, 
Physica E 40, 2642 (2008). 

A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. 
Novoselov, and A.K. Geim, Rev. Mod. Phys. 81, 10 (2009). 
S. Dubois, A. Lherbier, A. Botello-Mendez, J.-C. Charlier, 
in preparation. 

39 A. Lherbier, S.M.-M. Dubois, X. Declerck, S. Roche, Y.M. 
Niquet, and J.-C. Charlier, Phys. Rev. Lett. 106, 046803 
(2011). 

R. Haydock, V. Heine, and M.J. Kelly J. Phys. C : Solid 
State Phys. 8, 2591 (1975). 

The same broadening parameter (10 meV) has been used 
in both DOS calculations. 

T.O. Wehling, S. Yuan, A.L Lichtenstein, A.K. Geim, and 
M.L Katsnelson, Phys. Rev. Lett. 105, 056802 (2010). 



17 



A. Mesaros, S. Papanikolaou, C.F.J. Flipse, D. Sadri, and 
J. Zaanen, Phys. Rev. B 82, 205119 (2010). 

G. Trambly de Laissardiere, and D. Mayou, Modern Phys. 
Lett. B 25, 1019 (2011). 

A. Lherbier, S. Roche, O.A. Restrepo, Y.M. Niquet, A. 
Decorte, and J.-C. Charher, submitted for pubhcation 
(2011). 

V. Kapko, D.A. Drabold, and M.F. Thorpe, Phys. Status 
Sohdi B 247, 1197 (2010). 

E. Holmstrom, J. Fransson, O. Eriksson, R. Lizarraga, B. 
Sanyal, S. Bhandary, and M.I. Katsnelson, Phys. Rev. B 
84, 205414 (2011). 

S. Roche, and D. Mayou, Phys. Rev. Lett. 79, 2518 (1997). 
D. Mayou, and S. Khanna, J. Phys. I France 5, 1199 
(1995). 

^° D. Mayou, Europhys. Lett. 6, 549 (2007). 

F. Triozon, J. Vidal, R. Mosseri, and D. Mayou, Phys. Rev. 
B 65, 220202 (2002). 

S. Roche, Phys. Rev. B 59, 2284 (1999). 

H. Ishii, F. Triozon, N. Kobayashi, K. Hirose, S. Roche, 
C.R. Physique 10, 283 (2009). 

P.A. Lee, and T.V. Ramakrishnan, Rev. Mod. Phys. 57, 
287 (1985). 

P.A. Lee, and D.S. Fisher, Phys. Rev. Lett. 47, 882 (1981). 
F. Evers and A.D. Mirhn, Rev. Mod. Phys. 80, 1355 
(2008). 

F. Triozon, S. Roche, and D. Mayou, RIKEN Rev. 29, 73 
(2000). 

This also explain why the use of the same value of Dma^ 
has been preferred to normalize the two equations in order 
to clearly see this spurious effect. 

N.H. Shon, and T. Ando, J. Phys. Soc. Jpn. 67, 2421 
(1998). 

T. Ando, International Journal of Modern Physics B 21, 
1113 (2007). 

S. Roche, B. Biel, A. Cresti, F. Triozon, Physica E (in 
press), DOI: 10.1016/j-physe.2011.06008 (2011). 



T. Stauber, N.M.R. Peres, and F. Guinea, Phys. Rev. B 
76, 205423 (2007). 

J. Yan, and M.S. Fuhrer, Phys. Rev. Lett. 107, 206601 
(2011). 

J.H. Chen, C. Jang, S. Xia, M. Ishigami, and M.S. Fuhrer, 
Nature Nanotech. 3, 206 (2008). 

S.V. Morozov et al. , Phys. Rev. Lett. 100, 016602 (2008). 
E. Akkermans, and G. Montambaux, Mesoscopic Physics 
of Electrons and Photons, Cambridge University Press, 
ISBN 978-0521855129 (2007). 

J. Bang, and K.J. Chang, Phys. Rev. B 81, 193412 (2010). 
Y.V. Skrypnyk, and V.M. Loktev, Phys. Rev. B 83, 085421 
(2011). 

G. Schubert, and H. Fehske. [a?Xlv: 1109. 6439 vl (2011). 
R. Grassi, T. Low, M. Lundstrom, Nano Lett. 11, 4574 
(2011). 

D.C. Elias, R.R. Nair, T.M.G. Mohiuddin, S.V. Morozov, 
P. Blake, M.P Halsall, A.C. Ferrari, D.W. Boukhvalov, 
M.I. Katsnelson, A.K. Geim, and K.S. Novoselov, Science 
323, 610 (2009). 

'^^ A. Bostwick, J.L. McChesney, K.V. Emtsev, T. Seyller, 
K. Horn, S.D. Kevan, and E. Rotenberg, Phys. Rev. Lett. 
103, 056404 (2009). 

'^^ N. Leconte, J. Moser, P. Ordejon, H. Tao, A. Lherbier, A. 
Bachtold, F. Alsina, CM. Sotomayor Torres, J.-C. Char- 
lier, and S. Roche, ACS NANO 4, 4033 (2010). 
J. Moser, H. Tao, S. Roche, F. Alsina, CM. Sotomayor 
Torres, and A. Bachtold, Phys. Rev. B 81, 205445 (2010). 
N. Leconte, A. Lherbier, F. Varchon, P. Ordejon, S. Roche, 
and J.-C. Charlier, Phys Rev. B 84, 235420 (2011). 
A. Uppstu, K. Saloriutta, A. Harju, M. Puska, and A.-P. 
Jauho. larXivnrriI70r4 8vl (2011), submitted to Phys. Rev. 
Lett. 

J. Haskins, A. Kinaci, C. Sevik, H. Sevincli, G. Cuniberti, 

and T. gagin, ACS NANO 5, 3779 (2011). 

S. Wu, and F. Liu, arXiv:1001.2057^1 (2010). 

W. Li, Y. He, L. Wang, G. Ding, Z-.Q. Zhang, R. W. Lortz, 

P. Sheng, and N. Wang, Phys. Rev. B 84, 045431 (2011). 



